function gesol = PriceRealBonds_Model20220926(parms,gesol_in,FLAG_CALC_ENTROPY)

    if nargin<3
        FLAG_CALC_ENTROPY = true;
    end

    gesol = gesol_in;
    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    
    StateTransitionProbs_Nxzz = reshape(permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]),...
        [NN*Nx*Nz, Nz]); % ((N,L,z),z') format
    SDF_Nxzz = reshape(gesol.spd1period,[NN*Nx*Nz, Nz]); % ((N,L,z),z') format
    log_SDF_Nxzz = log(SDF_Nxzz);
    P_next = zeros(NN*Nx*Nz,Nz); % S'=S'((N,L,z),z')
    E_dlogM_next = zeros(NN*Nx*Nz,Nz);
    
    % law of motion
    x_next_Nxzz = gen_law_of_motion_x(parms);
    Nnext_Nxz = gen_law_of_motion_N(parms,gesol);
    
    [N_meshgrid,x_meshgrid] = meshgrid(parms.grid_N(:),parms.grid_x(:)); % needed for interp2
    
    Ts = (0:60)'; NT = numel(Ts);
    Ps = zeros(NN,Nx,Nz,NT);
    
    % initialize
    Ps(:,:,:,1) = 1;
    
    if FLAG_CALC_ENTROPY
        E_dlogM = zeros(NN,Nx,Nz,NT); % E_t[log(M(t+T)) - log(M(t))]
        E_dlogM(:,:,:,1) = 0;
    end
    
    for n = 2:NT
        
        % interpolate: transpose data input to account for difference
        % between meshgrid and ndgrid
        for iz_next = 1:Nz
            P_next(:,iz_next) = interp2(N_meshgrid,x_meshgrid,Ps(:,:,iz_next,n-1)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iz_next));
            if FLAG_CALC_ENTROPY
                E_dlogM_next(:,iz_next) = interp2(N_meshgrid,x_meshgrid,E_dlogM(:,:,iz_next,n-1)',...
                    Nnext_Nxz(:),x_next_Nxzz(:,iz_next));
            end
        end
        Ps(:,:,:,n) = reshape(sum(StateTransitionProbs_Nxzz.*SDF_Nxzz.*P_next,2),...
            [NN Nx Nz]);
        if FLAG_CALC_ENTROPY
            E_dlogM(:,:,:,n) = reshape(sum(StateTransitionProbs_Nxzz.*(log_SDF_Nxzz + E_dlogM_next),2),...
                [NN Nx Nz]);
        end
        
    end
    
    % store + convert to (N,H,z,T) format
    gesol.real_bonds.format = 'P=P(N,lambda,z,maturity)';
    gesol.real_bonds.maturities = Ts(2:end);
    gesol.real_bonds.prices = Ps(:,:,:,2:end);
    
    if FLAG_CALC_ENTROPY
        gesol.entropy.horizon = Ts(2:end);
        gesol.entropy.E_dlogM = E_dlogM(:,:,:,2:end);
        gesol.entropy.CondEntropy = log(gesol.real_bonds.prices) - gesol.entropy.E_dlogM;
        gesol.entropy.note = 'CondEntropy = CondEntropy(N,lambda,z,horizon)';
    end

end